#!/usr/bin/python
# -*- coding: utf-8 -*-

# Copyright (c) 2011
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
#
# Author: Jesus Carrero <j.o.carrero@gmail.com>
# Mountain View, CA

from scipy import r_, array, zeros, ones

class FemBasisBase(object):

    __slots__ = ['m_b_type', 'm_mass', 'm_deg_of_freedom', 'm_support']

    def __init__(self, b_type, _dof, support):
        """   """

        self.m_b_type = b_type
        self.m_deg_of_freedom = _dof
        self.m_support = support
        self.m_mass = None
        self._mass()

    def _mass(self):
        """ virtual method to call mass matrix."""

        assert False

    def dof(self):
        """ degrees of freedom   """

        return self.m_deg_of_freedom

    def shape_func_support(self):
        """ return extrems defineing domain for shape functions. """

        return self.m_support

    def get_mass_mat(self):
        """ get mass matrix.  """

        return self.m_mass

    def integrate(self, i, j):
        """ integrate two shape functions in reference interval."""

        return self.m_mass[i, j]


class LagrangeP1(FemBasisBase):

    """ compute the inner product
      ------- xi \in [0, 1]
  """

    def __init__(self):
        """ Instantiate an obj of type P1.  """

        FemBasisBase.__init__(self, 'lagrangeP1', 2, array([0, 1]))

    def l1(self, points):
        """  L_1(pointsi) = 1.0-pointsi """

        return 1. - points.reshape(points.size, 1)

    def l2(self, points):
        """ L_2(pointsi) = pointsi """

        return points.reshape(points.size, 1)

    def _mass(self):
        """ compute mass matrix. """

        self.m_mass = zeros((2, 2), 'float64')
        self.m_mass[0, :] = array([1., 1 / 2.])
        self.m_mass[1, :] = array([1 / 2., 1.])
        self.m_mass /= 3.

    @classmethod
    def get_grad_grad(cls):
        """ return grad-grad matrix.  """

        grgr = ones((2, 2), 'float')
        grgr[0, :] = array([1., -1.])
        grgr[1, :] = array([-1., 1.])
        return grgr


    @classmethod
    def get_grad_disp(cls):
        """ returns grad-disp matrix.  """

        grds = ones((2, 2), 'float')
        grds[0, :] = array([-1., 1.])
        grds[1, :] = array([-1., 1.])
        return grds / 2.

    def values(self, points):
        """ evalue shape function at quad points.
            quad formulas are defined in the interval [-1, 1].
        """
        qpoints1 = (points + 1.) / 2.
        basis_at_points = [self.l1(qpoints1)]
        basis_at_points.append(self.l2(qpoints1))
        basis_at_points = array(basis_at_points)
        basis_at_points = basis_at_points.reshape((2, points.size))
        return basis_at_points

    @classmethod
    def derivatives(cls, points):
        """ evaluate derivatives at quad points."""

        return r_[-ones((1, points.size)), ones((1, points.size))]

    @classmethod
    def second_der(cls, points):
        """ second derivative at quad points."""

        return zeros((2, points.size))

class LagrangeP2(FemBasisBase):

    """ compute the inner product
      ------- xi \in [-1, 1]
  """

    def __init__(self):
        """ Instantiate an object to P2 type. """

        FemBasisBase.__init__(self, 'LagrangeP2', 3, array([-1, 1]))

    @classmethod
    def l1(cls, points):
        """ L_1(pointsi) = -0.5*pointsi*(1.0-pointsi) """

        val = -.5 * points * (1. - points)
        return val.reshape(val.size, 1)

    @classmethod
    def l2(cls, points):
        """  N_2(pointsi) = 1.0-pointsi**2 """

        val = 1. - points ** 2
        return val.reshape(val.size, 1)

    @classmethod
    def l3(cls, points):
        """ N_3(pointsi) =  0.5*pointsi*(1+pointsi) """

        val = .5 * points * (1 + points)
        return val.reshape(val.size, 1)

    def _mass(self):
        """ compute mass matrix.  """

        self.m_mass = zeros((3, 3), 'float64')
        self.m_mass[0, :] = array([4., 2., -1.])
        self.m_mass[1, :] = array([2., 16., 2.])
        self.m_mass[2, :] = array([-1., 2., 4.])
        self.m_mass /= 30.

    @classmethod
    def get_grad_grad(cls):
        """ compute grad-grad at reference element.  """

        grgr = ones((3, 3), 'float')
        grgr[0, :] = array([7., -8., 1.])
        grgr[1, :] = array([-8., 16., -8.])
        grgr[2, :] = array([1., -8., 7.])
        return grgr / 3.

    @classmethod
    def get_grad_disp(cls):
        """ compute matrix grad-displacment.  """

        grds = ones((3, 3), 'float')
        int_ni = array([1., 4., 1.]) / 3.
        intx_ni = array([-1., 0., 1.]) / 3.

        grds[0, :] = -.5 * (int_ni - 2 * intx_ni)
        grds[1, :] = -2. * intx_ni
        grds[2, :] = .5 * (int_ni + 2 * intx_ni)
        return grds.T

    def values(self, points):
        """ evaluate shape function at points."""

        basis_at_points = zeros((3, points.size))
        basis_at_points[0, : ] =  self.l1(points).flatten()
        basis_at_points[1, : ] =  self.l2(points).flatten()
        basis_at_points[2, : ] =  self.l3(points).flatten()
        return basis_at_points

    @classmethod
    def derivatives(cls, points):
        """ evaluate at the derivatives of the shape fucntions. """

        tmp = points.reshape((1, points.size))
        return r_[r_[tmp - .5, -2. * tmp], tmp + .5]

    @classmethod
    def second_der(cls, points):
        """ second derivative of shape functions. """

        tmp = ones( (3, points.size) )
        tmp[1, :] *= -2.
        return tmp


